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ABSTRACT 

We investigate the large-scale clustering and gravitational interaction of baryons and 
dark matter (DM) over cosmic time using a set of collisionless N-body simulations. Both 
components, baryons and DM, are evolved from distinct primordial density and veloc- 
ity power spectra as predicted by early-universe physics. We first demonstrate that such 
two-component simulations require an unconventional match between force and mass res- 
olution (i.e. force softening on at least the mean particle separation scale). Otherwise, the 
growth on any scale is not correctly recovered because of a spurious coupling between the 
two species at the smallest scales. With these simulations, we then demonstrate how the 
primordial differences in the clustering of baryons and DM are progressively diminished 
over time. In particular, we explicitly show how the BAO signature is damped in the 
spatial distribution of baryons and imprinted in that of DM. This is a rapid process, yet 
it is still not fully completed at low redshifts. On large scales, the overall shape of the 
correlation function of baryons and DM differs by ~ 2% at z = 9 and by 0.2% at z = 0. 
The differences in the amplitude of the BAO peak are approximately a factor of 5 larger: 
10% at z = 9 and 1% at z = 0. These discrepancies are, however, smaller than effects 
expected to be introduced by galaxy formation physics in both the shape of the power 
spectrum and in the BAO peak, and are thus unlikely to be detected given the precision of 
the next generation of galaxy surveys. Hence, our results validate the standard practice of 
modelling the observed galaxy distribution using predictions for the total mass clustering 
in the Universe. 

Key words: cosmology:theory - large-scale structure of Universe. 



1 INTRODUCTION 



JLJ .In the standard paradigm of cosmological structure forma- 
• _ 'tion, primordial density perturbations are a result of quantum 
fluctuations amplified by cosmic inflation. At these very early 
<-h ■times, baryons and dark matter (DM) density fields have the 
_^_'same phases and amplitudes - each of them has a fluctua- 
tion spectrum following a power law with an index close to 
unity. However, subsequent interaction with the radiation field 
breaks the initial similarity, creating a scale-dependent growth 
that is different for baryons and for DM. 

On scales smaller than the horizon and prior to recom- 
bination, baryons couple to photons through Compton scat- 
tering. Radiation pressure opposes gravity and inhibits the 
growth of density perturbations. The balance is not perfect 
however, thus generating oscillations on sub-horizon scales 
in the density, temperature and size of perturbations in the 
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baryon-photon fluid. On even smaller scales, free-streaming 
and an imperfect coupling between baryons and photons pro- 
gressively damps the amplitude of these oscillations. In con- 
trast to baryons, DM particles do not directly interact with 
photons, and are thus mainly affected by gravity. DM den- 
sity fluctuations can grow freely, and are only halted by the 
Meszaros effect on scales smaller than the horizon at the 
matter-radiation equality. The physics describing these inter- 
actions is understood at high precision and is able to describe 
at very high accuracy the patterns of temperature fluctuations 
obs erved in the cosmic m icrowave background radiation (see 
e.g. iHu fc Dodelsonll2002l . for a review). 

After recombination, baryons decouple from the photons 
leading to a drop in sound speed by ~ 5 orders of magnitude 
and an associated drop in the Jeans mass of ~ 14 orders of 
magnitude. From now on, the evolution of perturbations in 
baryons and DM is dominated by the same physics and the 
growth is almost entirely determined by gravity until much 
later times (when hydrodynamical interactions become im- 
portant). The initial conditions for the two components after 
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recombination are, however, very different. The coupling be- 
tween baryons and photons has prevented the baryons from 
falling into the density fluctuations present in the DM fluid. 
The power spectrum of density fluctuations in baryons and 
DM are thus genuinely distinct. In particular, baryonic acous- 
tic oscillations (BAO) dominate the baryon density power 
spectrum but they are almost non-existent in the DM dis- 
tribution. 

At later times, the gravitational coupling between DM 
and baryons will reduce such differences. The DM distribu- 
tion gradually obtains a BAO signal, while the amplitude of 
BAOs in the baryons clustering is reduced. This process is 
commonly assumed to be finished at low redshift, i.e. baryons 
and DM are assumed to have identical spatial distributions 
(equal to that of the total mass field) on large scales. A nat- 
ural corollary of this is that the BAO should be detectable 
in the galaxy distribution. This has indeed been achieved ob- 
servationally with inc r easing accuracy (e.g. ICole et al.ll2005l : 
lEisenstein et all 120051 : iBlake et all 1201 li ; iBeutler et al*ll201lD 
and further measurements have been proposed using virtu- 
ally any known tracer of the matter density fi eld in the Un i- 
verse: including the galaxy distribution (e.g. ICooravl I2002T ). 
gala xy clusters (e.g Angulo et al.l 2005), the Ly — a forest 
(e.g. IWhite et al)l20ld : lKitaura et al.ll2012f) or 21 cm emissio n 
background from galaxies at low redshifts ( Chang et alj|2008h . 
from the epoch of re-ionisation at z ~ 10 (IMao fc Wul 120081; 
Rhook et all 120091 ). and even using Supernovae tZha n et al.1 



2008). 



However, the details of the process in which the DM 
(and thus the total matter) distribution acquires the BAO 
signature still remain relatively unexplored. At high redshift 
or large scales, the interaction between baryons and DM 
particles can be followed accurately by perturb ation theory 
(jSomogvi fc Smithll2010l : lBernardeau et al.ll2012T l. However, to 
explore low redshifts and small scales in the mass field, N- 
body simulations are essential, since they provide the most 
accurate and faithf ul predictions in the nonlinear regime (see 
iKuhlen et aTll2012l , for a recent review). Unfortunately, to our 
knowledge, no N-body simulation has been performed to ad- 
dress this topic. This is an important issue, since a precise 
understanding of the BAO signal in the z < 10 Universe is 
required to interpret accurately the high precision measure- 
ments that will be carried out over the next decade. 

In this paper, we directly follow the gravitational inter- 
action of DM particles and baryons, from z = 130 up to the 
present day. For this, we perform N-body simulations of two 
interacting fluids with different primordial density and veloc- 
ity fluctuations: one representing the DM field, and another 
representing the baryons. We provide details of these simula- 
tions and the numerical set-up in §2. With these simulations 
in hand, in §3, we explore the evolution of the large-scale clus- 
tering of baryons and DM, with particular emphasis on the 
evolution of the BAO peak. We discuss our results and con- 
clude in 64. 



2 N-BODY SIMULATIONS 
2.1 Initial conditions 

We generate the initial position and velocity f or our simula- 
tion p articles at z = 130 using the Music code (|Hahn fc Abell 
1201 ll ). and adopting a set of cosmological parameters con- 
sistent with the published measurements of the WMAP7 



data release l|Komatsu et al.lfeoill ). Explicitly: f2 m = 0.276, 
Q.A = 0.724, Q b = 0.045, h = 0.703, a s = 0.811 and spectral 
index n a = 0.961. 

We compute the primordial power spectra for baryons 
a nd DM using a linear Boltz mann solver code similar to that 
of iMa fc Bertschingerl (|l9951 ) where residual baryon-radiation 
interaction effects become small. The velocity field is irrota- 
tional and thus fully described by the velocity divergence, so 
that the initial conditions are fully specified by the power spec- 
tra of the overdensity 8 and the velocity divergence 6: Ps c , 
Pe c , Ps B , Pe B , where the subscript 'C stands for CDM and 
'B' for baryons. 

These power spectra are used together with the Zel'dovich 
app roximation (i.e. f irst order Lagrangian Perturbation The- 
ory [ZeTDovich| [l97(J) in Music to generate the initial condi- 
tions. The gravitational potential, whose gradient appears in 
both the particles displacement and velocity, is herein replaced 
by four potentials generated by the respective density power 
spectra and velocity divergence spectra. Hence, we set e.g. for 
the initial baryon positions and velocities 



x s = q + V<I>fl, v B = V* fl , 
where q is the Lagrangian coordinate and 

$s(k) oc g(k) k- 2 ^Ps B {k,t), and 



tfs(k) oc g(k)k~\ Pe B (k,t) 



(1) 

(2) 
(3) 



are the respective potentials and Q is a real-valued Gaussian 
random field of zero mean and unit variance. We refrain here 
from applying 2LPT as we start at a very high redshift. In 
addition, we note that in this way the streaming velocity be- 
tween baryons and DM is included self-consistently, but not 
its non-linear impact onto the spectra until z = 130, which 
however is negligibly small on the large scales we consider. 

We also note that at a redshift of 130 there remains a small 
contribution of the radiation energy density to the Friedmann 
equation. We ignore this in our simulations and linear theory 
calculations, setting fl r — at z ^ 130, but not before. 

The baryonic particles are placed on a staggered initial 
mesh with respect to the DM particles. The required phase 
shift of the noise field is computed in Fourier space. Staggering 
is necessary to minimise a spuriously tight cou pling between 
the two particle types (cf. Yoshi da et al.l I2003L but also our 
discussion in Section |2.3[) . We also note that the initial con- 
ditions for all our simulations use the same Gaussian white 
noise field, and thus, have identical random phases allowing a 
more direct comparison between them. 



2.2 Specifics of the N-body simulations 

We use two sets of cosmological N-body simulations to study 
the gravitational coupling of baryons and DM. In the first 
set, all matter is represented by a single fluid sampling the 
total matter power spectrum. The second set contains simu- 
lations with two distinct fluids, one representing baryons and 
the other DM, both of which have different primordial den- 
sity and velocity fluctuations, as predicted by early universe 
physics. Each set consists of two simulations with different 
box sizes: i) L = 1 /i -1 Gpc, sufficient to measure reliably the 
BAO signal; and ii) L — 250ft _1 Mpc, with which we focus on 
smaller scales and on the coupling between baryons and DM 
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Figure 1. A test of the effect of force resolution in N-body sim- 
ulations. The plot compares the ratio of the power spectrum of 
baryons and CDM at z = 19 for four runs with different force reso- 
lution. Solid and dotted lines show results for two different softening 
lengths corresponding to 100 /i — x kpc and 52 h~ J kpc (1/10-th and 
1/20-th of the mean inter-particle separation, respectively). Dashed 
and dot-dashed lines show results for two runs where forces are com- 
puted using a particle-mesh (PM) method, where the size of each 
grid cell is ~ 4h _1 Mpc and ~ lh _1 Mpc (corresponding to 1 and 
0.25 times the mean inter-particle separation). The red solid line 
displays the expectation of linear perturbation theory, whereas the 
blue solid line shows the results of a simulation where the forces are 
adaptively smoothed. 

in the nonlinear regime. In all cases, the DM and baryonic 
field was each represented with 1024 3 particles, which for the 
1 h~ 1 Gpc box implies a DM particle mass of 5.97x 10 10 /i _1 M 
and a baryon particle mass of 1.16 x W 10 h' 1 Mq. For the 
smaller simulation box these values are 64 times smaller: 
9.3 x 10 s /i _1 M Q and 1.81 x 10 8 h _1 M Q , respectively. 

We follow the non-linear evolution of the two fluids using 
a memory-efficient v ersion of the P-Gadget3 Tree-PM code 
(|Springel et al.ll2005l ) described in lAngulo et all (|2012D . In all 
cases, we use a mesh to compute forces of dimension equal to 
the particle mesh so to avoid spurious features at the beating 
frequency of the two meshes (|Angulo et al.l [20081). Note that 
we have neglected all hydrodynamical interactions and thus 
baryons behave as a collisionless fluid which is a good approx- 
imation on the large scales and for the processes we explore 
in this work. 



2.3 Force softening and discreteness effects 

The choice of the force resolution with which the N-body cal- 
culations are carried out is a sensitive issue for our work. Figfl] 
compares the relative difference in the clustering of baryons 
and DM, after they have been evolved from z — 130 up 
to z = 19, for different numerical configurations and force 
resolutions. Due to computational reasons, for this test we 
employ simulations with the same mass resolution as our 
L = 1000/i _1 Mpc run, but on a 8 times smaller volume 
L — 500/i -1 Mpc. We observe that the results we obtain from 
the different runs vary widely. Solid (diamonds) and dotted 
(triangles) lines show results for a Plummer softening length 



equal to 1/20-th and 1/10-th of the mean inter-particle sepa- 
ration (the typical softening scales adopted in state-of-the-art 
cosmological N-body simulations). Despite being a standard 
set-up, they overestimate the strength of the coupling between 
our two particle species. The same occurs for our 3rd test case, 
denoted by a dot-dashed line (crosses) , where we have disabled 
the short-range interactions but increased the size of the PM 
mesh to be equal to 1/4-th of the mean inter-particle separa- 
tion of each particle type. As in the previous test cases, the ra- 
tio of the power spectra still appears to be significantly higher 
than the prediction of linear theory (indicated by the solid red 
line). Theses calculations depart even further at lower redshift. 
A similar b ehaviour of incorrect l arge-s cale growth has been 
reported bv lO'Learv fc McQuinnl l|2012l ). 

We can track this issue to the collissionality seen in the 
runs with high force resolution, which stems from the discreti- 
sation of the density fields. On small scales, the force gen- 
erated by N-body particles only approximately represents an 
homogeneous and continuous force field. This fact causes the 
particles to quickly couple on small scales, which affects the 
evolution even on large scales, increasing the rate at which 
differences in the clustering of baryons and DM dissipate. 

We also note that we did not observe a significant impact 
of the choice of initial particle distribution on large-s c ale be - 
haviour. While lYoshida et al.l (|2003l ) and iNaoz etall (|201ll ). 
e.g., have argued for two randomly displaced glass distribu- 
tions for the two species, we find that any choice of 'staggered' 
initial distribution - where separations between particles from 
the two species are locally maximized - reduces, but does not 
avoid, the discreteness effects. Regardless of the initial parti- 
cle distributions, a very large force softening h as to be chosen 
to obt ain t he correct growth. T he reason why lYoshida et al.l 
(2003) (and lHahn fc Abelll201ll ) do not observe a similar spu- 
rious coupling is most likely due to the significantly smaller 
boxes investigated there, where strong non-linear growth dom- 
inates quickly over the numerical artefacts. 

It is important to note that the reason why discreteness 
effects are not as evident in common o ne- or two-component 
particle simulations as in our case (e.g. Irlamana et al.ll2002l ) 
is because these simulations typically start from identical per- 
turbation spectra, so that the spurious coupling cannot be 
easily diagnosed using the power spectrum. We argue, how- 
ever, that in two-component simulations (such as SPH+DM 
simulations), it should still appear as an additional binding 
energy between the two fluids that has to be overcome by pres- 
sure forces at late times. The impact is even less clear though, 
and harder to quantify, in standard one-fluid CDM calcula- 
tions, but we expect it to appear as an artifical population of 
small-scale of low- mass halos, somewhat similar to those seen 
in warm DM simulations (e.g. lWang fc Whitdl2007h . Whether 
this has any sizeable and unforeseeable consequence for the 
results of simulations remains to be investigated. 

We see in Fig. [T] that the only set-up that correctly re- 
produces the linear growth is the dashed line (squares), where 
the force has been smoothed on scales below the mean inter- 
particle separation for each component separately (or 1.25 
times the inter-particle separation for both components to- 
gether). (Despite this, the total mass power spectrum always 
displays the correct evolution.) Therefore, we will adopt this 
numerical configuration for the simulations used in the remain- 
der of the paper, i.e. using a grid of 1024 3 points to solve the 
Poisson equation. 

We note that the use of an adaptive softening length (e.g. 
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Ilannuzzi fc Dolagl 1201 lh can also alleviate the artefacts de- 
scribed above. The red line in Fig. fT] shows the results of a 
run where the forces between baryon and DM particles (not 
however the intra-species forces) are softened adaptively using 
an SPH kernel with a width set by the distance to the 32nd 
neighbour (see Apendix A of lSpringel et al.ll200ll . for a discus- 
sion of this standard feature of Gadget). This, however, comes 
at the expense of suppressing the growth of small scale pertur- 
bations even more than in our low force resolution PM runs. 
Note that adaptive softening is also inherent in all c osmolog- 
ical a daptive me sh refinement codes such as Ramses (|T evssicr 
l2002h and ENZO (|Brvan fc Normanlll997l : lO'Shea et al.ll2004h . 
so that a spurious coupling between the two species is not ex- 
pected to occur. For these codes, the numerical evolution of 
the baryons on a grid comes however at the expense of nu- 
merical diffusion at the grid scale, leading to a comparable 
suppression of high -fc modes in the baryons (see Figure 23 of 
lHahn fc Ab"ell201ll ) as those in the adaptive SPH run. 



3 THE CLUSTERING OF BARYONS AND 
DARK MATTER 

Having identified a robust numerical setup, in this section, 
we present predictions for the large-scale clustering of baryons 
and dark matter as measured in N-body simulations. We start 
showing results in Fourier space (§3.1), then move to configu- 
ration space where we center our discussion on the BAO peak 
(§3.2). 



3.1 Power Spectrum 

We begin by presenting in Fig. [5] the power spectra, P(k), of 
DM and baryons in real space, measured from our simulations 
at different output redshifts. Measurements were performed 
by mapping the particle distribution, using a clouds-in-cell 
scheme, onto a 1024 3 grid and then Fast Fourier Transforming 
this field. We correct the effects of the assignment scheme by 
dividing each mode by the Fourier transform of a cubical top- 
hat, but we do not subtract a Poisson shot-noise term. Here we 
show results from our large- and small-box simulations. The 
minimum wavenumber plotted is 27r/1000 = 0.0062 /iMpc" 1 
and 27r/250 = 0.025 /iMpc" 1 for each of them. 

The top panel of Fig.[2]compares the growth of total mass 
power spectrum from our runs with that predicted by linear 
perturbation theory, which is displayed as solid black lines. We 
note that we have evolving the z = 130 linear theory power 
spectrum without considering the interaction of photons and 
baryons (effectively setting Q, r — 0), which allows a more di- 
rect comparison with our simulations, where only gravitational 
interactions are considered. 

We can see that our simulations closely reproduce the 
expected linear growth, differing only at the 4% level: The 
growth of the fundamental mode, from z = 130 until z = 0, is 
10438.5, whereas the linear theory prediction is 10015.6. This 
close agreement supports the correctness of our numerical cal- 
culation. On small scales, the expected nonlinear growth dom- 
inates, and the measured z — power spectrum is a factor of 
~ 20 larger than linear theory expectations at k = 1 ftMpc -1 . 

In the two bottom panels of Fig. [2] we plot the ratio be- 
tween the power spectra of baryons and DM. This highlights 
the differences in the overall shape as well as in the BAO sig- 
nature, visible in the range 0.1 < (fc/ftMpc" 1 ) < 0.3. At all 



redshifts, the curves are systematically different from unity, 
which implies that the overall shape of the power spectrum of 
baryons and DM is different, even on very large scales. Den- 
sity perturbations in the baryon density field are smaller than 
those in the DM at all times, resulting from the extra suppres- 
sion produced by radiation pressure before recombination. 

At the starting redshift of our simulations [z = 130), the 
amplitude of the power spectrum of baryons is 40% of that of 
the DM component, even on scales as large as the turn-over 
(k — 0.02 /iMpc - ). On smaller scales, the difference is con- 
stant down to the smallest scales that our simulations resolve. 
At later times, the gravitational interaction between baryons 
and DM particles couples the perturbation fields at all scales. 
As a result, both fields are homogenised and initial differences 
are progressively reduced. At z — 35 the baryon to DM P(k) 
ratio is ~ 0.85%, and ~ 0.95% at z = 9. At the latter redshift, 
fluctuations have been amplified by a factor of 100, but the 
differences between the baryon and DM power spectra are still 
approximately scale independent. 

Gravitational interaction continues at lower redshifts, 
but, due to the effect of dark energy, perturbations grow more 
slowly than at higher redshifts. By z — the baryons power 
spectra show a 1% suppression compared to that of DM, al- 
most the same difference present at z = 2. The ratio now 
displays a scale dependence - small scales approach unity at a 
slower rate - which can be interpreted as the DM field having 
experienced more nonlinear evolution and mode coupling than 
the baryonic field. This could be a residual effect of the DM 
power spectrum having higher amplitude than that of baryons. 

The amplitude of the oscillatory behaviour seen in the 
lower panels of Fig[2] arises from differences in the BAO fea- 
ture. The bigger their amplitude, the more dissimilar the BAO 
are in the two fluids. At the starting redshift, oscillations are 
large, i.e. BAOs are very weak in the DM but very notorious 
in the baryons. (For baryons, the amplitude of BAO is ~ 30% 
of that of the power spectrum, but for DM, it is only ~ 5%.) 
At lower redshifts, gravity will damp the BAO in the baryons 
and increase their amplitude in the DM. However, as in the 
case of the overall shape of the power spectra, even at z = 0, 
there are still some residual differences. In the next subsection 
we will investigate the evolution of the shape and amplitude 
of BAO signal in more detail. 

Fig. O shows a detailed comparison of our results with 
the prediction of linear perturbation theory. The top/bottom 
panel shows the ratio between the power spectrum of 
DM/baryons to that of the total mass, divided by the same 
ratio predicted in linear theory. Thus, departures from unity 
indicate where linear theory is not able to predict the rela- 
tive differences, with respect to the total mass, present in the 
spatial distribution of DM or baryons. Firstly, we note that de- 
partures in this ratio are much smaller, and appear on smaller 
scales, than absolute deviations of each component with re- 
spect to linear theory. Here, differences are below 0.5% for the 
baryons and below 0.1% for the DM, to be compared with one 
order of magnitude deviations in their individual amplitude 
(c.f. Fig. [2]). Hence, linear theory predicts this quantity at a 
remarkable accuracy. We note that our findings quali t ativel y 
agree with the analytical results of So mogyi fc Smith! l|201Ch . 
who, using renormalised perturbation theory extended to a 
multi-fluid case, also reported a positive deviation from unity 
on small scales for DM, and negative deviation for baryons. 

Overall, from this figure we get a picture complementary 
to our earlier discussion: Nonlinear evolution produces a faster 
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Figure 2. The growth of density fluctuations in the baryon and in dark matter fields, as measured by the spherically averaged power spectrum 
in real space. We show results for 13 redshifts, as indicated by the annotations in the figure, and for two different simulations containing the 
same number of particles but on boxes of side 1000 h _1 Mpc or 250 h _1 Mpc. Their predictions overlap over the range k = [0.03 — 0.3]/iMpc _1 . 
In the top panel, coloured lines show the measured total mass power spectrum in our simulations whereas dashed, dotted and solid lines 
show predictions of linear perturbation theory for the baryon, CDM and total mass power spectra, respectively. The lower panels show the 
ratio between the measured power spectra from baryons to that from CDM. Note that this ratio is not affected by the noise that arises from 
the finite sampling of large modes present in the simulation box. No Poisson shot-noise has been subtracted. Note the suppression of baryon 
perturbations relative to CDM perturbations at large k that is due to non-linear effects. 



growth in DM perturbations relative to that in baryons. For 
this reason, the amplitude of the baryon spectrum is overes- 
timated, whereas that of DM is underpredicted. This results 
in an overprediction of the baryon-to-DM power spectra ratio, 
which increases at smaller scales. Linear theory predicts only 



a 1% difference between the clustering of DM and baryons by 
z — 0, whereas in our simulations we measure a 50% larger dis- 
crepancy. These departures are small, but nevertheless exem- 
plify that gravitational evolution can be followed at some de- 
gree by linear or higher-order perturbation theories, but once 
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Figure 3. Comparison between simulation results and linear per- 
turbation theory predictions for the relative difference between the 
power spectrum of baryons and CDM to the total mass field (top 
and bottom panel, respectively). We show results from our lh~ Gpc 
box for k < 0.3 hMpc -1 , and from our 250/i -1 Mpc for smaller 
scales. Coloured lines show results at different rcdshifts, as indi- 
cated by the legend. 



the density contrast reaches values well above unity, numerical 
simulations are needed to properly and accurately follow the 
growth of structure. 



3.2 The BAO peak in the Correlation Function 

We now complement the Fourier-space results by presenting 
the correlation function, of baryons, DM and that of 

the total mass, as measured in our simulations. We compute 
correlation functions in Fourier space using a mesh of 1024 
points per dimension, which provides an accurate estimation 
of correlation functions our > 20/i~ 1 Mpc for our l/i _1 Gpc 
box. This approach is considerably faster than a direct pair 
count when computing correlation functions on large scales in 
catalogues containing a large number of particles. 

Fig. 3] shows the correlation functions computed in this 
way for both baryons and DM in our 1 ft -1 Gpc simulation. The 
top panel displays the measurements in real space while the 
bottom panel shows the results in redshift space using a plane- 
parallel approximation (i.e. the effect of peculiar motions in 
distance estimators is included by considering an observer at 
infinity). We show results on scales r — [60 — 130]/i _1 Mpc to 
focus on the BAO peak (appearing at ~ 110/i~ 1 Mpc). Each 
measurement has been divided by the square of the growth 
facto r, and additionally by the Kaiser "boost factor" (|Kaiserl 
1987), in the case of redshift-space results. Therefore, if the 
growth of a given component is completely linear and scale- 
independent, then curves at different redshifts would coincide 
in each panel. We see that this is indeed a good first-order ap- 
proximation, though systematic differences exist. We discuss 
this next. 

In the top-left panel of Fig. [4] we can see how the BAO 



peak is gradually imprinted in the DM distribution. At z = 
130, the correlation function is still close to a power-law with 
just a relatively small bump at 110/i _1 Mpc. However, in an 
extremely rapid process, the BAO emerges and reaches almost 
its full amplitude by z ~ 30, only ~ 100 Myrs after. During 
this period, we also see a small scale-dependent suppression of 
the correlation function at scales r < 80 /i _1 Mpc, to accommo- 
date for the growing peak. The subsequent evolution happens 
at a lower rate, and the maximum relative amplitude of the 
BAO is reached at z ~ 6, the moment in which mild nonlinear 
coupling of independent Fourier modes start to have a notice- 
able effect, decreasing aga in the amplitude of BAO peak and 
broadening its shape (e.g. lAngulo et aHl2005l ). 

The history of baryons, displayed in the middle panels, is 
the opposite. At the starting redshift the correlation function 
is dominated by the BAO peak, it is ~ 4 times larger than the 
broad-band shape of the correlation function, and a factor of 
~ 2 larger than the maximum relative amplitude ever present 
in the total mass field. This, however, quickly changes, with 
the amplitude of the peak decreasing linearly with redshift. As 
in for the DM case, density fluctuations in the baryons also 
show a scale-dependent growth, but in the opposite direction. 

The evolution in baryons and DM exactly compensate 
each other to give an almost perfectly scale-independent 
growth of fluctuations in the total mass density, as can be 
seen in the right panels of Fig. 0] This behaviour is only bro- 
ken at low redshift by non- linear mode couplin g , which smears 
out t he BAO peak (e.g. lAngulo et al.l 120081 : ISanchez et al] 
2008, and references therein). The evolution is compensated 
since perturbations do not grow independently, but are linked 
through a common gravitational potential that determines an 
identical acceleration field for both components. Thus, gravity 
acts as a homogenizer of the fields. In fact, although by z — 
there are residual differences, if we let the simulation run into 
the future, then eventually both fields will be indistinguish- 
able from each other on large scales. (We note that the same 
effect could be achieved by a larger value for as.) 

The bottom panel of Fig.[4]shows the spherically averaged 
clustering in redshift space (a quantity closer to that observed 
by spectroscopic galaxy surveys) for baryons, DM and for the 
total mass. For both components, and at all redshifts, we can 
see that the measured correlation functions are not simply a 
scaled version of their real-space counterpart, as one might 
naively expect in linear theory. 

At high redshift, the redshift-space enhancement in the 
overall correlation function is smaller than linear theory ex- 
pectations for DM, but larger for baryons. At low redshifts, 
predictions are closer to our measurements, and the respec- 
tive curves therefore overlap better in this plot. At the same 
time, the BAO peaks get damped as a consequence of nonlin- 
ear contributions to peculiar velocities. At high redshifts, this 
is also true for baryons but not for the DM field, which sees 
an increase in the amplitude of the BAO peak. 

All the differences described above are summarised and 
quantified in Fig. [5J which displays the correlation function 
for baryons and DM relative to that of the mass for differ- 
ent redshifts and for both real and redshift space, averaged 
over two scales. In the top panel, averaged between 60 and 
75/i _1 Mpc, capturing changes in the overall amplitude. In 
the bottom panel, between 105 and 110fo -1 Mpc, so it shows 
changes in the amplitude of the BAO peak. 

We see that the correlation function amplitude at z = 10, 
relative to that of the total mass, is ~ 1% higher for the case of 
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Figure 4. Real-space (top) and redshift-space (bottom) correlation functions in our simulations at different redshifts, focusing in the BAO 
peak. Measurements are shown for baryons, dark matter and total mass (panels (a), (b) and (c)), and for linear bins of Ar = 3h — x Mpc. 
Results were scaled by the square of the appropriate growth factor, and by the Kaiser boost factor in the case of the redshift-space correlation 
functions. 



DM, and ~ 3% smaller for baryons. These figures are reduced 
to a sub-percent level by z — 0. The differences are largest 
for the BAO peak in the baryons: ~ 10% excess at z = 10, 
which is reduced to about 1 percent by z = 0. This effect 
is smaller than those introduced by galaxy formation physics 
(Angulo et al. 2013, in prep.), which could cause up to ~ 10% 
effects. 

Redshift space distortions (RSD) reduce the amplitude of 
the BAO peak in the baryons at all redshfts, relative to that 
in real space. For DM, they enhance the BAO peak at high 
redshifts, but damp it after z ~ 1. The observed behaviour 
is explained by the fact that the gravitational potential that 
originates the coherent velocity flows is generated by the total 
mass distribution, which is not identical to either the baryonic 
or DM material but to their mass weighted average. Thus, the 
density-velocity relation expected in linear theory does not 
hold separately for either of the components. At high redshift 
this homogenises the field further at the BAO scale, whereas 
at low redshift nonlinear RSD dominate and diffuse the peak. 

3.3 Comparison with 1-fluid simulation 

Finally, we test the widespread assumption that the nonlinear 
evolution of perturbations in the density field of baryons and 
DM can be approximated as a single fluid representing the 



total mass field. For this, we run our two N-body simulations 
with exactly the same numerical configuration, but this time 
baryons and DM have identical primordial density and velocity 
power spectra (equal to that of the total mass), so that they 
effectively behave as a single fluid. With this we can test to 
what degree the nonlinear evolution of two interacting fluids 
is equivalent to that of a single fluid with the mass-weighted 
average power spectrum. 

Fig. [()] presents our results. In the top panel we show the 
ratio of the power spectra of our two- and one-fluid simula- 
tions. On large scales, both simulations are virtually indis- 
tinguishable, it is only in the nonlinear regime where differ- 
ences appear. The one fluid case underestimates the amount 
of nonlinear clustering in a roughly redshift- independent man- 
ner. The differences are small: at k < 2ftMpc _1 are less than 
0.1%, although they increase exponentially with the wavenum- 
ber k to reach percent level by k ~ 10. Unfortunately, our 
simulations lack the spatial resolution to pin down this more 
accurately. Our results agree with the analytical findings of 
ISomogvi fc Smith! {2010), who found a discrepancy of less than 
~ 0.3% at k = 1 /iMpc^ 1 between the mass power spectra from 
one- and two-fluid calculations at z = 0. 

In the bottom panel of Fig. [()] we display an analogous 
comparison but focusing on the BAO peak. Contrary to the 
Fourier-space behaviour, we see that the amount of nonlin- 
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Figure 5. The evolution of the average correlation function for 
baryons and DM relative to that of the total mass. In the top panel, 
we show the ratio of the correlation functions averaged between 60 
and 75 h~ *Mpc separations. Thus, this displays differences in the 
overall amplitude of the respective correlation function. In the bot- 
tom panel, we perform the average between 105 and 110 /i — 1 Mpc. 
Thus this can be regarded as the evolution in the amplitude of the 
BAO peak for baryons and DM. 
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Figure 6. Comparison of the power spectrum (top panel) and cor- 
relation function (bottom panel) predicted by our two-fluid simula- 
tions with those predicted by a single-fluid simulation, in which the 
initial conditions of particles were drawn following the total mass 
power spectrum at the starting redshift of our simulations. Lines of 
different colours display results at different redshifts, following the 
same convention as in Fig. [3] 



ear diffusion of the BAO peak is larger in the one fluid case. 
Nevertheless, the differences are extremely small: below 0.01%. 
These results validate the use of one-fluid simulations to study 
the large-scale distribution of mass, and hence of the BAO sig- 
nal, in the Universe. 



4 SUMMARY AND CONCLUSIONS 

We have performed a direct calculation of the gravitational 
coupling of baryons and dark matter, from z = 130 to z = 0, 
using N-body simulations. We focused primarily on large 
scales to explore how the BAO signature, initially present pri- 
marily only in the spatial distribution of baryons, is imprinted 
in the DM distribution. 

At the starting redshift, the density and velocity power 
spectra imprinted in baryons and DM differ considerably. At 
latter times, gravity acts as a diffusive agent of the primordial 
clustering differences. We find that the bulk of the differences 
are dissipated quickly, but also that this process is not fully 
completed by the present day. Initial differences in the power 
spectrum amplitude of ~ 40% at z = 130 are brought down to 
sub-percent level by z — 0. Interestingly, the BAO peak does 
not evolve as quickly and differences are still above 1% at low 
redshifts. 

On large scales, our results are consistent with linear 
perturbation theory: the ratio between the power spectra of 
baryons and DM decreases in a roughly scale-independent 
manner at high redshifts. On smaller scales of k > 0.2 /iMpc -1 , 
nonlinear evolution breaks this and we predict a growth of the 
ratio between the baryon and the DM power spectrum that 
is slower than in linear theory. The size of these departures is 
small (below 0.5%), thus, in practice, linear theory is a very 
accurate predictor of all these effects. 

Additionally, we confirmed that the nonlinear evolution of 
the total mass power spectrum can be accurately predicted by 
numerical simulations treating all mass as a single fluid with 
a single power spectrum. Most of these res ults are in qualita- 
tive a greement with the analytical work of lSomogyi fc Smith! 
i|2010T ). 

We conclude that assuming that baryons and dark matter 
have the same spatial distribution on the BAO scale is a good 
approximation given the accuracy with which it is expected 
to be measured by upcoming large galaxy surveys. However, 
we emphasize that the quality of this assumption is redshift- 
dependent, being worse at high redshifts and better at low 
redshifts. The amplitude in the BAO peaks differs at the 10% 
level at z = 10, but only at 1% at z = 0. The differences at low 
redshifts are smaller than the effects expected to be introduced 
by galaxy formation physics, but they could perhaps eventu- 
ally be detected in the future in, for instance, HI clustering 
from SKA measurements, or in another. If this is the case, our 
work indicates that at high redshift the bias of tracers that 
are sensitive to baryons instead of the total matter distribu- 
tion might include also a non-negligible baryon bias, though 
more work in this respect is needed, since the size of the ef- 
fect likely depends on the type of tracer and object selection 
criteria employed. 

Finally, we highlighted that force resolution is a critical 
issue to obtain accurate results for simulations with two flu- 
ids with distinct primordial density fluctuations. Discreteness 
effects in standard numerical configurations resulted in an arti- 
ficial coupling between baryons and DM particles, which prop- 
agated to large scales and lead to an artificially fast growth 
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of large-scale modes. We tracked this issue to the loss of col- 
lisionallity in the fluids, which, thus, could be solved only by 
softening gravitational forces below the mean inter-particle 
separation. The numerical artefacts arising from a high force 
resolution were evident in our runs, but they should also be 
present in standard N-body calculations. However, while the 
net effect is not clear, it clearly warns for a careful assess- 
ment of the robustness of current numerical simulations in 
the mildly non-linear regime. 
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